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ABSTRACT 

The physical relationship between the far-infrared and radio fluxes of star 
forming galaxies has yet to be definitively determined. The favored interpre- 
tation, the "calorimeter model," requires that supernova generated cosmic ray 
(CR) electrons cool rapidly via synchrotron radiation. However, this cooling 
should steepen their radio spectra beyond what is observed, and so enhanced 
ionization losses at low energies from high gas densities are also required. Fur- 
ther, evaluating the minimum energy magnetic field strength with the traditional 
scaling of the synchrotron flux may underestimate the true value in massive star- 
bursts if their magnetic energy density is comparable to the hydrostatic pressure 
of their disks. Gamma-ray spectra of starburst galaxies, combined with radio 
data, provide a less ambiguous estimate of these physical properties in starburst 
nuclei. While the radio flux is most sensitive to the magnetic field, the GeV 
gamma-ray spectrum normalization depends primarily on gas density. To this 
end, spectra above 100 MeV were constructed for two nearby starburst galax- 
ies, NGC 253 and M82, using Fermi data. Their nuclear radio and far-infrared 
spectra from the literature are compared to new models of the steady-state CR 
distributions expected from starburst galaxies. Models with high magnetic fields, 
favoring galaxy calorimetry, are overall better fits to the observations. These so- 
lutions also imply relatively high densities and CR ionization rates, consistent 
with molecular cloud studies. 

Subject headings: galaxies: individual (NGC 253, M82) — galaxies: starburst — 
gamma-rays: galaxies 
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Introduction 



The successful operation of the Fermi Gamma-ray Space Telescope provides our first 
opportunity to infer the potential origi ns of the dif fuse high energy emission from normal 
galaxies beyond the Milky Way (e. g., iButt 1 12010| ). Recent detections of diffuse gamma- 



ray emissio n from Local Group, starburst, and Seyfert 2 g alaxies demonstrate the system's 
sensitivity (lAbdo et al. I l201ll . l2010al : iLenain et al. 1 12010| ) . With further integration time 
beyond these detection studies, we can construct gamma-ray spectra in order to constrain 
models of the cosmic ray (CR) interactions that give rise to the emission in these sources. 

For starburst galaxies in particular, the normalization of the gamma-ray spectrum helps 
to constrain, among other properties, the volume-averaged gas density affecting the CR dis- 
tributions. With this and the radio normalization, the strength of the so-called "calorimeter" 
model of CR cooling may be determined for these sources. At issue is the discrepancy in 
massive starbursts between the minimum energy magnetic field strength, estimated from the 
radio flux, and the equip artition field strength, es timated from gas surface density and assum- 
ing hydrostatic balance (IThompson et al. II2006I ). Given th eir large surface densities and the 
linearity of the FIR-radio correlation over decades in flux, [Thompson et al. I (120061 ) argues 
for relatively high field strengths in starbursts. These values can be an order of magnitude 
above the minimum energy field estimates, and thus alter significantly our understanding of 
the energy budgets in starburst galaxies. 

In this work, the gamma-ray, radio, and FIR spectra of the archetypal starburst galaxies 
NGC 253 and M82 are quantitatively compared to self-consistent models of the various energy 
loss and emission mechanisms for CRs in normal galaxies. A broad parameter space is 
explored, including density and magnetic field strength, to constrain the physical conditions 
of these starburst nuclei that best agree with their full spectral energy distributions. 



Observations and Results 



The data were obtained using the Fermi Large Area Telescope (LAT) (lAtwood et al 



20091 ) . The LAT is a pair-conversion detector sensitive to photons with ener gies between 



rough ly 20 MeV and 300 GeV. The publicly available Fermi Science Tools (lAbdo et al. 



20091 ) were used, and the threads reviewed on the Fermi Science Center web sitqj for the 



binned likelihood analysis were followed. 



^http : //f ermi . gsf c . nasa . gov/ ssc/ data] 
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Almost 3.5 yr of Pass 7 LAT data were downloaded from 10° regions of interest around 
NGC 253 and M82. To avoid gamma-rays from the Earth's limb, a zenith cut of 100° was 
invoked in the data filtering step. The rock angle was constrained to be less than 52° and 
various other data flags and filters were set as suggested by the LAT team. With the large 
field of view of the LAT, the background must be modeled in order to ensure photons only 
associated with the target source are counted. The 2FGL catalog was used as an initial 
background, and then gttsmap was used to search for unmodeled sources that showed up 
with high significance. Two additional sources were identified within 5° of NGC 253, and 
three were found near M82. These new sources were incorporated into the background model 
with the likelihood routine gtlike. Thirty logarithmically-spaced energy bins per were used 
for the counts cube and binned exposure map. The spatial binning for both was set to 
0?1 per pixel. An additional 10° was added to the exposure map radius to account for the 
larger PSF at low energies. The P7SOURCE_V6 instrument response functions were used. 
The final energy binning was settled upon to maximize the signal-to-noise and detection 
significance, based on the test statistic (TS), in each bin. 

It should be noted that no source with a significant TS was detected above 100 MeV at 
the location of M81. However, the 0?61 separation between M82 and M81 is smaller than 
or comparable to the point spread function of the LAT below 1 GeV, so they are currently 
difficult to resolve. For this work, it was assumed that none of the flux from the position of 
M82 originates from M81. The FIR and radio luminosities of M82 are more than an order of 
magnitude higher than those of M81, so the gamma-ray contribution from M81 is expected 
to be negligible. 



For NGC 253, the integrated flux between 100 MeV and 100 GeV was (1.19±0.15) x 10 



-8 



ph cm"^ s"\ with a photon index of -2.24 ± 0.06, and TS of 113. For M82, a total flux 



of (1.44 ± 0.22) X 10^® ph cm~^ s^^ was found, with a photon in dex of — 2.25 ± 0.08 , and 
TS of 180. These results are consistent wi th initial detections (lAbdo et al. I l2010bl ) and 
the se cond LAT catalog (I Ab do et al. 1120111). They also extend to meet the TeV observa- 



tions (lAcero et al. II2009I : lAcciari et al. 1120091 ). Lightcurves show no detectable variability 



in eight- week bins. The differential spectra of both galaxies are shown in Figures [T] and [21 
and are summarized in Table [H A linear fit yields values that indicate the spectra are 
inconsistent with a simple power law. In particular, there may be an apparent turnover in 
the spectra below a GeV indicative of the expected CR proton interactions in these sources 
(IBlom. Paglione fc Carramiiianalll999l ). The following section presents models of gamma-ray 
and radio emission from the expected CR populations in starburst nuclei. 
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3. Modeling the Starburst Emission 



The likely origin of gamma-rays from starburst galaxies is the interaction of supernova- 
accelerated CRs with a large mass of dense, pri marily molecular, gas . The high temperatures 
of dense molecular clouds in starbursts (e. g., iPaglione et al. 1120041: ISakamoto et al. 1201 if) 



indicate significant penetration and ionization by CRs (ISuchkov. Allen. &: Heckman 



1993|). 



The CR electrons themselves emit synchrotron radiation, so radio spectra of starburst re- 
gions serve as a significant constraint as mentioned before. The imbedded, dusty star forming 
regions should radiate strongly in the FIR, and if the CR electrons lose most of their en- 
ergy within them, this would be a natural explanation for the FIR-radio correlation and 
calorimetry. However, the gamma-ray spectra are critical for helping to reduce degeneracies 
and ambiguities in the physical solutions. 

To quantify the starburst contribution to the gamma-ray and radio emission from star- 
burst galaxies, the steady-state CR distributions were ni odeled. The models are significantly 
updated ver sions of earlier studies ( IPaglione et al. Ill996l ). similar in appro ach to more rec ent 



work 



.e. g. 



Torres 1120041 ). The steady-state CR diffusion-loss equation ( jLongair Ill98ll ) is 



now integrated over energy, where the original work solved a time integration, requiring a 
slightly different Green's function. For electrons, as before, energy losses due to ionization, 
non -thermal bremsstrahlung, synchrotro n radiation, IC scattering, and escape are accounted 
for ( Paglione et al. I 19961: Torres 2004 , appendix). The full Klein-Nishina cross section is 



now included ( Schlickeiser 2002[ p. 88) for IC scattering since relativistic effects become 



important at the high energies to which Fermi and ground-based arrays are now sensitive. 
Bremsstrahlung losses are now included in the continuous CR energy loss rate rather than 
approximating it with an energy-independent CR loss timescale. Also the photon energy 



density for IC is now constrained using the starburst templates of ISiebenmorgen &: Kriigel 



(120071 ) fit to the IR spectra. Finally, positrons (see below) are now included. 



The CR protons are affected by ionization and escape losses, p lus pion produc t ion a bove 
a total energy of 1.22 GeV. The cross section parameterizations of iKamae et al. I (120061 ) are 
now employed for all pion species (vr^'"). Secondary CR electrons are created via knock-on 
or Coulomb interactions, and the decay of negatively charged pions. Secondary positrons 
are created through the decay of positive pions. The distributions of t hese secondary l epton s 



from charged pions are also determined using the parameterizations of iKamae et al. I ( 120061 ). 



Power law primary CR electron and proton injection spectra of the form Q{E) = KE^'^, 
in particles cm~'^ s~^ MeV~^, are assumed, where E is the particle t otal energy. T he previous 
model used a power law in kinetic energy, similar to other studies (ITorres 1120041 ). which has 
a conveniently simple overall normalization. The total energy is chosen in this work no t only 



to be consistent with other theoretical studies (ILacki. Thompson, fc QuataertI l2010l : iBell 
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19781 ). but as a deliberate contrast from the previous work to test the effect of a different 
injection spectrum. This choice mostly affects the proton spectrum at low energies, and thus 
the pion spectra, which, via secondary production, could alter both the radio and gamma-ray 
emission, particularly near the pion bump. 

This is strictly a one-zone model only simulating the expected nuclear emission from 
starburst galaxies. No diffusion to, or emission from, the extended disk is considered. There 
are a variety or reasons for this choice besides basic simplicity. Given that these sources 
are so near the flux and angular resolution detection thresholds already, the much lower 
emissivitie s from the disk are i nsignificant to Fermi. We originally modeled the disk of 
NGC 253 ( iPaglione et al. Ill996[ ) and indeed found its gamma-ray emissivity to be negligible 
compared to th e central s tarburst, in the sam e way as its radio emission is much reduced 
( lUlvestadl l2000l ) . Recently [Torres et al.l ( 120121 ) estimated more realistically the steep decline 
(by an order of magnitude) of the CR proton density outside the inner starburst in NGC 
253. The diffuse gamma-ray and rad io emission from the Milky Way Galactic center also 
outshine the disk (jCrocker et al.ll201ll ). Our nearest neighbor,the LMC, one of the very few 
extragalactic sources with resolved gam ma-ray emis sion, similarly only shows emission from 
its most massive site of star formation (IButt 1120101 ). 



Using a power law injection spectrum in total energy also requires a new normalization. 
As usual, the normalization K is determined by the supernova rate, \E', starburst volume, 
V , and the efficiency of power transfer between the supernovae and CRs, ?7, such that 



g(E)EH„ dE = r]P^/V 



The supernova energy, P = 10^^ ergs. A minimum CR kinetic energy, E^in > 2mpv1 ~ a few 



MeV , is required for acceleration in a shock front given shock velocities of lO^kms ^ (IBell 



19781 ). and is much smaller than the proton rest mass energy. The relative normalization 
between the proton and el ectron CR distributions at high energies then de p ends on the 



inject ion spectral index s (IBell I Il978l : iLacki et al. I l2010l : iPersic et al.l l2008l : ISchlickeiser 



2OO2I . p. 472), and as such is no longer an independent free parameter: 



NJK = K/m,)(^-i)/2 



(2) 



Typical values of s = 2.0-2.4 yield proton-to-electron ratios of 50-200, which are consistent 
with Milky Way observations and effectively const rain the choi c e of s . Note that this is a 
sir nplification which differs from the treatment of iLacki et al. I (120101 ) (see also discussion 



m 



Persic et al.ll2008l ) in order to reduce the number of free parameters. This normalization 



mostly impacts the secondary contribution to the CR electron+positron flux, and introduces 
variations of at most factors of ~ 2 from other treatments. 
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A convection loss timescale Tc = 1 Myr is assumed, which is consistent with esti- 



mates of starburst w inds and CR transport time (or length) scales (iHeesen et al. 1 12009 



Murphy et al. Il2008[ ). So the total energy-dependent escape timescale is the combined dif- 



fusion and convection terms, or 



(3) 



In sum, the model variables are s, the volume-averaged gas (H2-I- H) density n, magnetic field 
strength B, diffusion time scale tq, and photon field energy density Uph (from the integrated 
starburst SED template). 



3.1. Model Results: Steady-State CR Distributions 

Figure [3] shows the results from a few models broadly indicating the effects of varying 
density and magnetic field strength on the steady-state CR electron+positron distributions. 
Higher densities increase the ionization losses at low energies, leading to an overall flatter, 
or more curved, spectrum. High magnetic fields increase synchrotron radiation losses, sub- 
sequently decreasing the CR electron distribution, particularly at high energies. As a result, 
high energy emission from inverse Compton (IC) scattering and non-thermal bremsstrahlung 
are reduced. Losses from IC scattering are less important; the overall results are quite in- 
sensitive to changes in Uph or the FIR spectrum template. Increased losses due to a shorter 
convection timescale (tc) are also relatively small, primarily lowering and slightly steepening 
the CR proton sp ectrum around TeV energies. Longer diffusion timescales (higher tq, see 



Torres et al.ll2010l ) result in slightly harder and higher CR spectra above 10 GeV or so. This 
effect is also comparatively small, reaching factors of just up to a few above a TeV, and 
negligibly changing the radio spectrum. 

Altering the injection slope s affects the steady-state spectral slopes accordingly, but 
can also substantially change the contribution of secondary particles because of its effect on 
the proton-to-electron ratio. For example, for flatter slopes, secondary particles contribute 
more to the lepton distribution (Eq. [2]), causing a more pronounced bump around 0.1 GeV 
in the spectrum. Note that this secondary bump is very subtle compared to the primary 
spectrum, so adjusting it by factors of a few changes the total lepton spectrum very little. 
Thus the results are not very sensitive to whether a power law injection spectrum in total 
or kinetic energy is chosen, nor to the specific treatment of Np/Nf,. 

The CR proton spectrum depends primarily on the density. As with the electrons, 
higher densities slightly flatten the proton spectrum at low energies from ionization losses. 
However, high energy losses are also increased via pion production, so the overall effect is 
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that the proton distribution decreases, and slightly flattens, with density. The effect of using 
an injection spectrum that is a power law in total, rather than kinetic, energy mostly appears 
in the resultant pion distributions (and thus the secondary electron+positron distributions). 
The pion spectrum is marginally lowered (by factors of ~ 2) at the 0.1 GeV bump, but is 
unaffected at high energies. This small difference seems to have a negligible effect on the 
resultant gamma-ray and radio emission spectra as the secondary contribution is slight, and 
the high energy pion distributions are unaltered. 



3.2. Model Results: Gamma-ray and Radio Spectra 

Radio and gamma-ray spectra were generated from the steady-state CR distributions. 
Gamma-rays are created by IC scattering and non-thermal bremsstralung from primary 
and secondary CR electrons-|-positrons, plus the decay of neutral pions from CR protons 
interacting with ambient gas. Radio spectra are generated by synchrotron radiation from 
these CR leptons, thermal free-free absorption and emission, plus warm dust emission at 
millimeter wavelengths. We neglect TeV absorption and pair production for now as their 



expec ted impact is relatively small and difficult to discern with current sensitivities ((Torres 



20041 ). 



The gamma-ray spectrum of a star forming galaxy above a GeV is dominated by pion 
decay from CR protons interacting with the interstellar medium (FigureHj). At lower gamma- 
ray energies, bremsstrahlung and/or IC scattering from CR electrons is important. The IC 
contribution is very sensitive to the high energy tail of the CR lepton distribution, and 
contributes most when there is a sufficient density of very high energy CRs, that is, when 
the magnetic field - and associated energy losses - is low. In fact, sensitive MeV data 
would be very useful to constrain the TeV CR lepton distributions in these galaxies. Higher 
densities enhance not only pion decay emission, but also non-thermal bremsstrahlung, since 
the secondary CR electron and positron distributions increase with density given that they 
too arise from CR protons (via charged pion decay). Therefore, the gamma-ray spectrum 
normalization is most sensitive to density, with some dependence on magnetic field strength 
below a GeV. 

The radio spectrum normalization is very sensitive to the magnetic field. Higher gas 
densities tend to reduce the radio signal somewhat by fiattenin g the CR electron distri 



bution. The fiatteniii g, or curvature, in the radio spectrum (jWilliams &: Bower 1 12010 



Thompson et al. II2006I ) is a result of ionization reactions affecting mostly lower energy CR 
electrons. Higher densities also enhance the bump feature in the secondary distributions, 
adding to the curvature. To a much lesser extent, the far- infrared (FIR) emission from dust 
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also affects the radio and gamma-ray spectra by providing a deep reservoir of seed pliotons 
for IC scattering. So a liigli FIR luminosity increases CR leptonic losses, consequently dimin- 
ishing slightly the non-thermal bremsstrahlung and synchrotron emission, while enhancing 
the IC gamma-ray flux. The FIR-radio correlation hints at the interrelatedness of all of these 
processes in star-forming galaxies. Neither the radio nor gamma-ray emission is very sensi- 
tive to the diffusion timescale. Together the radio and gamma-ray normalizations effectively 
constrain n and B. 



3.3. Model Fits to NGC 253 and M82 



The Spitzer IR spectra of NGC 253 and M82 are shown in F i gures [5] and [6] with 
some closely matching starburst SEDs from ISiebenmorgen fc Kriigel I ( 120071 ). The photon 
densities from these models are 171 and 274 eV cm~^ for NGC 253 and M82, respectively. 
The assumed supernova rates are 0.1 and 0.08 yr~^ for M82 and NGC 253, respectively, and 
their assumed distances are 3.9 and 3.5 Mpc. 

The radio and gamma-ray spectra were simultaneously fit using hundreds of model 
runs over a broad parameter space. G iven the possibility that the magnetic fields in such 
strong starbursts may be quite high ([Thompson et al. I l2006l ) . field strengths up to mC 
were included. Despite the small effect, we also investigated very short and long diffusion 
timescales, from 0.1 to 30 Myr since recent gamma-ray studies o f supernova-cloud i nteractions 
indicate that long timescales reproduce the observations best (ITorres et al.ll2010l ). The best 
fitting models were determined using the minimum x^- The radio data were selectively 
chosen to include only the central starburst disks (data from within the inner arcmin, or 
central 300-500 pc where the molecular and mm continuum emission are strongest), not 
the more extended winds/ halos. In starburst galaxie s , this extended emis s ion tends to have 



a mu ch steeper spectrum (jWilliams fc Bower 



2010; Heesen et al. 2009: Elmouttie et al 



19971 ) and is derived from a more evolved CR population rather than one directly interacting 
with the central starburst and contributing to the pion decay gamma-rays. The volumes for 
these galaxies were assumed to be disks of 150 pc radius and 70 pc thickness. This volume 
may underestimate their full molecular extent, especially in M82, but this distinction is not 
so critic al here; fiuxe s are generated from the modeled emissivities by multiplying by the 
volume (ITorres II2004I ). thus canceling with the 1/V term in the normalization (Equation [T]). 



The best fit models are shown with the data in Figures There is a large mismatch 
in both the noise and number of data points between the radio and gamma-rays. Thus the 
scatter in the gamma-ray models is understandable since the normalization is dictated mostly 
by the radio data. To stress the physical implications of the new Fermi data, fits were done 
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artificially increasing the weighting of the gamma-ray data by the ratio of the numbers of 
radio and gamma-ray points. Very similar solutions result, though over a narrower range of 
parameter space. 



Lacki et al.l ( 1201 ll ) compared the first Fermi catalog detections and TeV data to a 



collection of published model predictions and concluded that they were for the most part 
broadly consistent with each other (many of these models were constrained merely by non- 
detections at GeV energies and so considerably overestimated the observations). The scatter 
in the models in their Fig. 1 mirrors the 3-0" envelope of fits in Figs. 

Normalizing the radio and gamma-ray models to match the spectra yields power transfer 
efficiencies r] f» 6% for both galaxies. All fits within 3-a of the mini mum result in efficien- 
cies w ell below 20%, which is consistent with theoretical predictions ( iMarkiewicz. Drury fc Volk 
1990h . 



Solutions with low values (Figures [TO] and [TT]) do include low magnetic fields around 
the minimum energy value, B ~ 150/iG, but with very low gas densities of ~ 100 cm~^. 
Based on studies of dust and molecular gas in these galaxies, volurne-averaged densities of 
at least several hundre d H2-I-H cm~^ are expected (jPaglione et al. 1120041 : ISakamoto et al. 



201 ll : IWeiB et al.ll200ll ). Models with higher densities imply higher magnetic fields of several 
hundred /iG, approaching the cal orimeter mod e l exp ectations. Valid solutions also result 



for short confinement timescales. [Torres et al.l (120101 ) found that a diffusion coefficient of 
D = 10^^ cm^ s~^ is strongly constrained by Fermi and TeV observations of supernova- 
cloud interactions. Given tq ^ R'^/2D, and setting R to the starburst radius of 150 pc, 
yields tq = 8 Myr. We restrict the following discussion just to models where n > 300 cm~^ 
and Tq = 10 Myr. 



4. Discussion 



These result s generally favor the hig h magnetic field calorimeter model for starburst 



CR populations (IThompson et al. I l2006l ). The implication is that the CR cooling rate. 



primarily via synchrotron radiation and ionization losses, is faster that the escape rate. 
Figure [12] shows the loss timescales for electrons using models consistent with the radio 
and gamma-ray data for NGC 253. Here the CR lifetime is simply approxim ated as either 



the energy loss time scale E/b{E), where b{E) is the CR energy loss rate ([Torres [ [2004 



Paglione et al. [ [19961 ) . or the steady-state lifetime, N{E)/Q{E). Even for lower densities 



and magnetic fields, the electron loss timescales are shorter than the escape timescale. Even 
protons appear mostly calorimetric, especially for low energies and high densities. So it 
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would seem that the condition that tioss <^ Tesc is satisfied for electrons, particularly for 
higher values of n and B and given that diffusion timescales are likely longer than 1 Myr 
(recall that depends weakly on tq). Therefore the minimum energy estimate for magnetic 
field strength may indeed underestimate the true value in starburst galaxies. 

Another argument in favor of the high field (high density) solutions comes from a rough 
examination of the transport and interactions of CRs during their lifetimes within the star- 
burst region. CRs encounter a column of primarily molecular gas whi ch has been measured 



from CO studies to have beam-average d values of up to 10 H2 cm (iPaglione et al. 112004 



Paglione et al. II2OOII : IWeiB et al.ll200l[ ). where typical single-dish beam sizes encompass the 
entire central starburst region (300-600pc). One may approximate this column density as 
nAl, where Al is akin to the mean free path of the CR. This path length is then of order 
100 pc for low average gas densities (30 cm~^) and around a few pc for high values (10^ 
cm~^). The former result implies that the typical CR leaves the starburst region within its 
lifetime, rather than losing the majority of its energy within it, contradicting the steady- 
state assumption and calorimetry. But it is hard to understand independent studies of 
very warm cloud core temperatures in starbursts if CRs interact in such a limited way with 
molecular cloud cores (explored below). The high density (and high magnetic field) solu- 
tions imply that CRs lose the majority of their energy on spatial scales of order the sizes 
of dense molecular cloud clumps or cores. This result is much more consistent with model 
expectations and observations from photodisso ciation and/or cosmic ray dominated regions 
([Matsushita et al. II2OIOI: iPapadopoulos I I2OIOI ). the calorimeter model, the high measured 
cloud temperat ures (IPanuzzo et al.ll20101) . and the size scales of Galactic molecular clumps at 



such densities (IMcOuinn et al 



l2002h . e the measured volume 

densities and peak column densities of individual , ^ 30 pc scale cloud c omplexes in NGC 
253 and M82 (10^2-24 H2 cm'^ and lO^-^ cm'^) JSakamoto et al. II2OI1I : IWeiB et allboOlh . 



Comparable magnetic field strengths, and volume and column densities, have been inferred 
across the full ~ 40 p c extent of the massive star-forming cloud Sgr B2 in the Galactic center 



fICrutcher et al.lll996h . 



Further justification for the high density and magnetic field solutions may be derived 
from the expected ionization rates from penetrating CRs in starburst molecular clouds. Ion- 
ization rates, ( = Ce+Cpy ^-^e estimated using the best fit models for CR electron-|-positron and 
proton distributions, and the cross sections and procedures from lPadovani. Galli. fc Glassgold 
(120091 ). Only the H2 cross sections are used in this analysis since they are the most signifi- 
cant. The steady-state CR fiuxes are extrapolated down to the H2 ionization energy of 15.6 
eV for the integration. Note that this extrapolation overestimates the true ionization rate 
since other absorption p rocesses such as H2 exc itation and dissociation may occur as gas col- 
umn densities get large. iPadovani et al. I (120091 ) indicates that at low energies (< 1 MeV and 
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< 100 MeV for electrons and protons, respectively), the CR spectrum becomes rather flat, 
with a plateau that declines with the encountered column density. Deep within clouds, the 
CR proton contribution dominates over CR electron ionization given the apparently more 
gradual decline of (p with column density. 

This analysis yields median values of C ~ 6xl0~^^ s~^ for NGC 253 and M82 (Fig- 



ure [T3|l . This CR ionization rate is a few orders of magnitude above the e stimate in M82 
based on CR heating and the observed cloud temper atures, C ~ 4xl0~^^ ( ISuchkov et al. 



1993h. and factors of 1 -10 above the rate for ULIRGs ( Papadopoulos 2010). Taking ( from 



Suchkov et al. I ( 119931 ) as the fiducial rate necessary to achieve the observed cloud kinetic 



temperatures in M82 and NGC 253, and reducing t he mean C from this wo rk by that mag- 



nitude requires a column density of ~ 10 H2 cm ( iPadovani et al. Il2009[ ). consistent with 



the observed column density. Further, Figure [T3] indicates trends between ( and density and 
magnetic field strength. Solutions for lower n and B require more a ttenuation of the low 



energy CR flux to match the fiducial (. From iPadovani et al. I ( l2009l ). column densities of 



up to 10^^ H2 cm~^ or more are needed to reduce ( by the required factors of 10^ (Figured 
Such column densities exceed the gas mass of their entire starburst regions, which further 
supports the higher density solutions, and thus the high B calorimeter model. At very 
high densities however, lower columns are needed to reach the fiducial (, and Al = N-a^/n 
can yield path lengths much less than a pc, which is unreasonably small. Therefore, the 
ionization rate analysis disfavors the extremely high and low density (and magnetic field) 
solutions. 

Finally, theoretical and observed scaling relations between B, and A^H2 and n, complete 
the argument against very low magnetic fields in these starbursts. While the radio and 
gamma-ray solutions indicate a degeneracy between n and B (Figures [TO] and [TT]) . which 
also roughly follows the expected sc aling relation B oc y/n, the magnetic field should also 
increase proportionately with A^H2 (jCrutcherl Il999l ). This scaling is shown in Figure [TH 
Again, magnetic fields of several hundred fiG are favored. 



5. Conclusions 

Self-consistent models of the steady-state CR distributions in starburst galaxies match 
both the radio and gamma-ray spectra of NGC 253 and M82. The models solve the CR 
diffusion-loss equation given a power law spectrum of supernova-accelerated particles, ac- 
counting for energy and particle losses due to ionization, synchrotron radiation, IC scatter- 
ing, bremsstrahlung, secondary particle generation, and escape. The radio normalization 
depends mostly on the magnetic field, with some sensitivity to density in mostly the curva- 
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ture of the spectrum. The gamma-ray spectrum, now observable with sufficient sensitivity 
with Fermi and TeV detections, is mostly normalized by the density, with some dependence 
on the magnetic field below a GeV. Therefore, with both radio and gamma-ray data, the 
densities and magnetic fields in a starburst region may be independently and less ambigu- 
ously constrained than from radio data alone. In particular, the question of calorimetry 
(high magnetic fields) in massive starbursts may be addressed. 

In NGC 253 and M82, th e se mo del results favor the high magnetic field solutions an- 
ticipated by [Thompson et al. I (120061 ). The higher densities also implied are consistent with 
molecular gas studies. Further, the ionization rates, approximated from the steady-state 
CR distributions with high magnetic fields and densities, also match expectations. These 
results confirm that CR penetration and heating is an important contributor to the warm 
temperatures observed in starburst galaxies. These starburst nuclei may be populated to 
a significant filling factor with giant molecular clouds resembling Sgr B2 in magnetic field 
(~ 500;uG), column density (10^^ II2 cm~^), and density (10^ cm~'^). 
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Fig. 1. — Differential gamma-ray spectru m of NGC 253. Fi lled symbols are Fermi data 
(this work), the open symbol is from HESS (lAcero et al. II2009I ) assuming a photon index of 
-2.2 ±0.3. 
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Fig. 2. — Differential gamma-ray spectr um of M82. Filled s ymbols are Fermi data (this 
work), open symbols are from VERITAS ( lAcciari et al. 1 120091 ) . 
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Fig. 3. — Normalized total (secondary and primary) electron+positron steady-state distri- 
butions. The thick line shows model results for a nominal photon field and relatively high 
gas densities and magnetic field strengths. The dashed line shows the effect of lowering the 
magnetic field, and the dotted line shows the effect of a low density. Note in the high den- 
sity model the slightly more pronounced bump around 0.1 GeV from enhanced high energy 
secondaries. For all models here, s = 2.2 and Tq = 10 Myr. 
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Fig. 4. — Sample model results. (Left) Expected radio synchrotron emission for nominal 
values of density, magnetic field, and photon density (black). Here s = 2.2, r = 10 Myr, 
and 7] = 0.02. The blue line is for a lower magnetic field, and the red line is for a lower 
density. The dot-dash line indicates a typical thermal bremsstrahlung level and the onset of 
free-free absorption at low frequencies (here Tighz = 0.04). The thermal absorption has not 
been applied to the synchrotron spectra shown here. (Right) Expected gamma-ray emission 
due to neutral pion decay (solid lines), non-thermal bremsstrahlung (dotted lines), and IC 
scattering (dashed lines). Colors are as before. The green line indicates a lower photon 
density which is indistinguishable from the nominal case except for IC scattering. 
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Fig. 5. — IR spectrum of NGC 253. The SED models shown are for a 350 pc radius starburst 

region, a ratio of OB luminosity to total luminosity of 0.6, total IR luminosity of 

Ay = 72 mag, and a hot spot density of 10^ cm~^ (thin line). The thick line is for a 
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Fig. 6. — IR spectrum of M82. The SED model shown is for a 350 pc radius star- 
burst region, a ratio of OB luminosity to total luminosity of 0.9, total IR luminosity of 
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Fig. 7. — Radio (left) and gamma-ray (right) spectra from the central few hundred pc 
of NGC 253. Model solutions within 3cr of the minimum are shown in grey with the 
minimum {x^ = 2.2) in black. The best solution given the priors tq = 10 Myr and 
n > 300 cm~^ is indicated with the dashed line {B = 350/iG, n = 1 0^ cm~^, s = —2.4, 
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Fig. 8. — Same as Figure [TJ but after increasing the weighting of the gamma-ray data by the 
ratio of the numbers of radio and gamma-ray points (minimum = 4.9). The best solution 
for the priors tq = 10 Myr and n > 300 cm~^ is indicated with the dashed hne {B = 350/iG, 
n = 10^ cm-3, s = -2.4, Uph = 200 cm'^, r/ = 0.05, = 5.3.). 
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Fig. 9. — Same as Figure [TJ but for M82 (minimum = 3.5). The best solution for the 
priors tq = 10 Myr and n > 300 cm~^ is indicated with the dashed hne [B = 450/iG, 
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Fig. 10. — Density and magnetic field solutions for the best fitting models for NGC 253. 
Here s = 2.2, Upu = 200 eV cm-=^, and tq = 10 Myr. Contours indicate the 67%, 90%, and 
99% confidence levels around the minimum = 2.5. 




Fig. 11.— Model solutions for M82 as in Fig [TOl Here s = 2.2, Uph = 270 eV cm-^, and 
To = 10 Myr. 
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Fig. 12. — CR electron and proton loss and escape timescales. Electron energy loss timescale, 
E/b{E) (short dashed lines); electron steady-state lifetime, N{E)/Q{E) (solid lines); electron 
escape timescale, Eq. [3] (long dashed line). Here s = 2.2, Uph = 270 eV cm~^, and tq = 1 
Myr. Upper curves are for slow loss rates, n = lOOcm"^ and B = 150/iG; lower curves are for 
fast loss rates, n = 10^ cm~^ and B = TOOyuG. Proton loss {N/Q only) and escape timescales 
are denoted with thick lines. Again, the upper and lower curves are for n = 100 and 10^ 
cm~^, respectively. 



-29- 



10 



-9 



10 



-10 



+ 

Q. 



-11 



10 



-12 



10 



-13 



T| 1 1 1 — I I I I I I 1 1 1 — I I I I I I 1 1 1 — I I I I, 

o 



o 

o 
o 



10^ 



oO 

o» ° 



Jj I I I I I I I I I ■ ■ ■ ■ 



— 1 — I — I — r-| 

o 



— 1 1 1 1 — I — I — n— 



o 

O o 
o o 

o fiO 

O °o 

0^ 

0^ 



§ 



_i I I I I I 



10^ 10^ 
n (cm"^) 



100 1000 
B(iiG) 



Fig. 13. — Ionization rates for M82 (circles) and NGC 253 (diamonds). Note that ( here is 
expected to overestimate the actual ionization rate given the unmodeled attenuation of the 
low energy CR fluxes as they pass through molecular clouds. 
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Fig. 14. — Column densities needed to attenuate the ionization rates in M82 to the fiducial 
value of 4xl0~^ ^ s~^. Dashed h ues indicate the range of column densities estimated from CO 



mass ( iPanuzzo et al 



201^. The thick line in dicates the relat ion between column density and 



obser vations (WeiB et al.ll2001l ). and the dash-dot line is deduced from the central dynamical 



magnetic field in Galactic molecular clouds (ICrutcherlll999l ). Filled symbols denote column 
densities resulting in CR path lengths 1 pc < A/ = Nn^/n < 200 pc. 
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Table 1. Fermi LAT Differential Spectra of NGC 253 and M82 
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